# Replication files for 
# Masterson and Lehmann
# "Refugees, Mobilization, and Humanitarian Aid: Evidence from the Syrian Refugee Crisis in Lebanon"

# File 2: Appendix Analysis
# Requires File 1 be run before File 2

####
# Table of Contents -------------------------------------------------------

# 1 Movers Statistics 
## 1.1 How many HHs moved from one village to another?  
## 1.2 How many HHs moved to a different altitude level
## 1.3 Examine whether there were any differences among movers between experimental groups

# 2 Tables
## 2.1 Contamination of the Control Group 
## 2.2 Additional Aid Binary
### 2.2.1 Linear analysis without controls -- Additional aid binary
### 2.2.2  Linear analysis with controls -- Additional aid binary
##  2.3 quadratic analysis without controls -- Additional aid binary 
#### 	2.3.1 quadratic analysis with controls -- Additional aid binary
## 2.4  -- Additional aid value
### 2.4.1 Linear analysis without controls -- Additional aid value
### 2.4.2  Linear analysis with controls -- Additional aid value
## 2.5 quadratic analysis without controls -- Additional aid value 
### 2.5.1 quadratic analysis with controls -- Additional aid value

# 3. Output Aid Contamination Results Table 

# 4. Sharing Results

# 5. Additional Aid Binary 
## 5.1 Linear analysis without controls -- Additional aid binary
###	5.2  Linear analysis with controls -- Additional aid binary
### 5.3 quadratic analysis without controls -- Additional aid binary 
###	5.4 quadratic analysis with controls -- Additional aid binary
## 5.5  -- Additional aid value
### 5.5.1 Linear analysis without controls -- Additional aid value
###	5.5.2  Linear analysis with controls -- Additional aid value
## 5.6 quadratic analysis without controls -- Additional aid value 
### 5.6.1 quadratic analysis with controls -- Additional aid value

# 6. Output Aid Contamination Results Table 

# 7. McCrary Tests Section 


####
# 1 Movers Statistics -------------------------------------------------------


#  total_movers.txt
# 1.1 How many HHs moved from one village to another?  
sum(data[, mover])
write(sum(data[, mover]), file = "output/total_movers.txt")


#  altitude_movers.txt
# 1.2 How many HHs moved to a different altitude level
sum(data[, AltMover])
write(sum(data[, AltMover]), file = "output/altitude_movers.txt")


#  moversCell11.txt  
#  moversCell12.txt   
#  moversCell21.txt   
#  moversCell22.txt   
# 1.3 Examine whether there were any differences among movers between experimental groups

######   control 0 to control 1
cell11 <- sum(data[, moverCtoC])

##### treatment 0 to control 1
cell21 <- sum(data[, moverTtoC])

#####  control 0 to treatment 1
cell12 <- sum(data[, moverCtoT])

###### treatment 0 to treatment 1
cell22 <- sum(data[, moverTtoT])

write(cell11, file = "output/moversCell11.txt")
write(cell12, file = "output/moversCell12.txt")
write(cell21, file = "output/moversCell21.txt")
write(cell22, file = "output/moversCell22.txt")  





####
#   2 Tables

####
# 2.1 Contamination of the Control Group --------------------------------------
#  otheraid_results.tex 

OtherAidValueMat <- cbind(data$Q47_1_2_USD, data$Q47_2_2_USD, data$Q47_3_2_USD, data$Q47_4_2_USD, data$Q47_5_2_USD,
                          data$Q47_6_2_USD, data$Q47_7_2_USD)

AidValue <- apply(OtherAidValueMat, 1, sum)

coefs_otherAid <- c()
sesCR_otherAid <- c()
p_valuesCR_otherAid <- c()


coefs_otherAidValue <- c()
sesCR_otherAidValue <- c()
p_valuesCR_otherAidValue <- c()


# Additional Aid Binary ---------------------------------------------------

# 2.2 Additional Aid Binary

####  2.2.1 Linear analysis without controls -- Additional aid binary

LinearHolder <- lm(data$Q46  ~ outcomesData_long[, "Treatment"] + outcomesData_long[, "WINTELEV_oct13_rescaled"] +  I(outcomesData_long[, "WINTELEV_oct13_rescaled"] * outcomesData_long[, "Treatment"])  )
#print(colnames(outcomesData_long)[i])
#print(clx(fm = LinearHolder,  cluster = outcomesData_long$cas_code))
clx(fm = LinearHolder,  cluster = outcomesData_long$cas_code)[2, 1]
clx(fm = LinearHolder,  cluster = outcomesData_long$cas_code)[2, 2]
clx(fm = LinearHolder,  cluster = outcomesData_long$cas_code)[2, 4]

clx(fm = LinearHolder,  cluster = outcomesData_long$cas_code)

coefs_otherAid[1] <- clx(fm = LinearHolder,  cluster = outcomesData_long$cas_code)[2, 1]
sesCR_otherAid[1] <- clx(fm = LinearHolder,  cluster = outcomesData_long$cas_code)[2, 2]
p_valuesCR_otherAid[1] <- clx(fm = LinearHolder,  cluster = outcomesData_long$cas_code)[2, 4]

##### 	2.2.2  Linear analysis with controls -- Additional aid binary


LinearCovHolder <- lm(data$Q46  ~ outcomesData_long[, "Treatment"] + outcomesData_long[, "WINTELEV_oct13_rescaled"] +  I(outcomesData_long[, "WINTELEV_oct13_rescaled"] * outcomesData_long[, "Treatment"]) +
                        outcomesData_long$age0to4 +
                        outcomesData_long$age5to12 +
                        outcomesData_long$age13to17 +
                        outcomesData_long$age18to59 +
                        outcomesData_long$age60plus +
                        outcomesData_long$Oct2013Pop +
                        outcomesData_long$disabled +
                        outcomesData_long$origin_gov_0 +
                        outcomesData_long$origin_gov_4 +
                        outcomesData_long$origin_gov_6 +
                        outcomesData_long$origin_gov_10 +
                        outcomesData_long$origin_gov_11 +
                        outcomesData_long$educ_le_priinc +
                        outcomesData_long$educ_pricom +
                        outcomesData_long$educ_middle +
                        outcomesData_long$educ_sec_tech_uni +
                        outcomesData_long$age_head +
                        outcomesData_long$monthinleb +
                        outcomesData_long$base_relatives +
                        outcomesData_long$base_friends 
)
#print(colnames(outcomesData_long)[i])
#print(clx(fm = LinearCovHolder,  cluster = outcomesData_long$cas_code))
clx(fm = LinearCovHolder,  cluster = outcomesData_long$cas_code)[2, 1]
clx(fm = LinearCovHolder,  cluster = outcomesData_long$cas_code)[2, 2]


coefs_otherAid[2] <- clx(fm = LinearCovHolder,  cluster = outcomesData_long$cas_code)[2, 1]
sesCR_otherAid[2] <- clx(fm = LinearCovHolder,  cluster = outcomesData_long$cas_code)[2, 2]
p_valuesCR_otherAid[2] <- clx(fm = LinearCovHolder,  cluster = outcomesData_long$cas_code)[2, 4]



####  2.3 quadratic analysis without controls -- Additional aid binary 



QuadHolder <- lm(data$Q46  ~ outcomesData_long[, "Treatment"] + outcomesData_long[, "WINTELEV_oct13_rescaled"] +  I(outcomesData_long[, "WINTELEV_oct13_rescaled"] * outcomesData_long[, "Treatment"]) +  outcomesData_long[, "WINTELEV_oct13_rescaled"]^2  +  I(outcomesData_long[, "Treatment"] * outcomesData_long[, "WINTELEV_oct13_rescaled"]^2) )
#print(colnames(outcomesData_long)[i])
#print(clx(fm = QuadHolder,  cluster = outcomesData_long$cas_code))
clx(fm = QuadHolder,  cluster = outcomesData_long$cas_code)[2, 1]
clx(fm = QuadHolder,  cluster = outcomesData_long$cas_code)[2, 2]

clx(fm = QuadHolder,  cluster = outcomesData_long$cas_code)


coefs_otherAid[3] <- clx(fm = QuadHolder,  cluster = outcomesData_long$cas_code)[2, 1]
sesCR_otherAid[3] <- clx(fm = QuadHolder,  cluster = outcomesData_long$cas_code)[2, 2]
p_valuesCR_otherAid[3] <- clx(fm = QuadHolder,  cluster = outcomesData_long$cas_code)[2, 4]





##### 	2.3.1 quadratic analysis with controls -- Additional aid binary


QuadCovHolder <- lm(data$Q46  ~ outcomesData_long[, "Treatment"] + outcomesData_long[, "WINTELEV_oct13_rescaled"] +  I(outcomesData_long[, "WINTELEV_oct13_rescaled"] * outcomesData_long[, "Treatment"]) + outcomesData_long[, "WINTELEV_oct13_rescaled"]^2  +  I(outcomesData_long[, "Treatment"] * outcomesData_long[, "WINTELEV_oct13_rescaled"]^2) +
                      outcomesData_long$age0to4 +
                      outcomesData_long$age5to12 +
                      outcomesData_long$age13to17 +
                      outcomesData_long$age18to59 +
                      outcomesData_long$age60plus +
                      outcomesData_long$Oct2013Pop +
                      outcomesData_long$disabled +
                      outcomesData_long$origin_gov_0 +
                      outcomesData_long$origin_gov_4 +
                      outcomesData_long$origin_gov_6 +
                      outcomesData_long$origin_gov_10 +
                      outcomesData_long$origin_gov_11 +
                      outcomesData_long$educ_le_priinc +
                      outcomesData_long$educ_pricom +
                      outcomesData_long$educ_middle +
                      outcomesData_long$educ_sec_tech_uni +
                      outcomesData_long$age_head +
                      outcomesData_long$monthinleb +
                      outcomesData_long$base_relatives +
                      outcomesData_long$base_friends 
)
#print(colnames(outcomesData_long)[i])
#print(clx(fm = QuadCovHolder,  cluster = outcomesData_long$cas_code))
clx(fm = QuadCovHolder,  cluster = outcomesData_long$cas_code)[2, 1]
clx(fm = QuadCovHolder,  cluster = outcomesData_long$cas_code)[2, 2]


clx(fm = QuadCovHolder,  cluster = outcomesData_long$cas_code)


coefs_otherAid[4] <- clx(fm = QuadCovHolder,  cluster = outcomesData_long$cas_code)[2, 1]
sesCR_otherAid[4] <- clx(fm = QuadCovHolder,  cluster = outcomesData_long$cas_code)[2, 2]
p_valuesCR_otherAid[4] <- clx(fm = QuadCovHolder,  cluster = outcomesData_long$cas_code)[2, 4]







# Additional Aid Value ----------------------------------------------------



# 2.4  -- Additional aid value

####  2.4.1 Linear analysis without controls -- Additional aid value


LinearHolderAidValue <- lm(AidValue  ~ outcomesData_long[, "Treatment"] + outcomesData_long[, "WINTELEV_oct13_rescaled"] +  I(outcomesData_long[, "WINTELEV_oct13_rescaled"] * outcomesData_long[, "Treatment"])  )
#print(colnames(outcomesData_long)[i])
#print(clx(fm = LinearHolderAidValue,  cluster = outcomesData_long$cas_code))
clx(fm = LinearHolderAidValue,  cluster = outcomesData_long$cas_code)[2, 1]
clx(fm = LinearHolderAidValue,  cluster = outcomesData_long$cas_code)[2, 2]

clx(fm = LinearHolderAidValue,  cluster = outcomesData_long$cas_code)

coefs_otherAidValue[1] <- clx(fm = LinearHolderAidValue,  cluster = outcomesData_long$cas_code)[2, 1]
sesCR_otherAidValue[1] <- clx(fm = LinearHolderAidValue,  cluster = outcomesData_long$cas_code)[2, 2]
p_valuesCR_otherAidValue[1] <- clx(fm = LinearHolderAidValue,  cluster = outcomesData_long$cas_code)[2, 4]


##### 	2.4.2  Linear analysis with controls -- Additional aid value


LinearCovHolderAidValue <- lm(AidValue  ~ outcomesData_long[, "Treatment"] + outcomesData_long[, "WINTELEV_oct13_rescaled"] +  I(outcomesData_long[, "WINTELEV_oct13_rescaled"] * outcomesData_long[, "Treatment"]) +
                                outcomesData_long$age0to4 +
                                outcomesData_long$age5to12 +
                                outcomesData_long$age13to17 +
                                outcomesData_long$age18to59 +
                                outcomesData_long$age60plus +
                                outcomesData_long$Oct2013Pop +
                                outcomesData_long$disabled +
                                outcomesData_long$origin_gov_0 +
                                outcomesData_long$origin_gov_4 +
                                outcomesData_long$origin_gov_6 +
                                outcomesData_long$origin_gov_10 +
                                outcomesData_long$origin_gov_11 +
                                outcomesData_long$educ_le_priinc +
                                outcomesData_long$educ_pricom +
                                outcomesData_long$educ_middle +
                                outcomesData_long$educ_sec_tech_uni +
                                outcomesData_long$age_head +
                                outcomesData_long$monthinleb +
                                outcomesData_long$base_relatives +
                                outcomesData_long$base_friends 
)
#print(colnames(outcomesData_long)[i])
#print(clx(fm = LinearCovHolderAidValue,  cluster = outcomesData_long$cas_code))
clx(fm = LinearCovHolderAidValue,  cluster = outcomesData_long$cas_code)[2, 1]
clx(fm = LinearCovHolderAidValue,  cluster = outcomesData_long$cas_code)[2, 2]

clx(fm = LinearCovHolderAidValue,  cluster = outcomesData_long$cas_code)

coefs_otherAidValue[2] <- clx(fm = LinearCovHolderAidValue,  cluster = outcomesData_long$cas_code)[2, 1]
sesCR_otherAidValue[2] <- clx(fm = LinearCovHolderAidValue,  cluster = outcomesData_long$cas_code)[2, 2]
p_valuesCR_otherAidValue[2] <- clx(fm = LinearCovHolderAidValue,  cluster = outcomesData_long$cas_code)[2, 4]


####  2.5 quadratic analysis without controls -- Additional aid value 


QuadHolderAidValue <- lm(AidValue  ~ outcomesData_long[, "Treatment"] + outcomesData_long[, "WINTELEV_oct13_rescaled"] +  I(outcomesData_long[, "WINTELEV_oct13_rescaled"] * outcomesData_long[, "Treatment"]) +  outcomesData_long[, "WINTELEV_oct13_rescaled"]^2  +  I(outcomesData_long[, "Treatment"] * outcomesData_long[, "WINTELEV_oct13_rescaled"]^2) )
#print(colnames(outcomesData_long)[i])
#print(clx(fm = QuadHolderAidValue,  cluster = outcomesData_long$cas_code))
clx(fm = QuadHolderAidValue,  cluster = outcomesData_long$cas_code)[2, 1]
clx(fm = QuadHolderAidValue,  cluster = outcomesData_long$cas_code)[2, 2]

clx(fm = QuadHolderAidValue,  cluster = outcomesData_long$cas_code)



coefs_otherAidValue[3] <- clx(fm = QuadHolderAidValue,  cluster = outcomesData_long$cas_code)[2, 1]
sesCR_otherAidValue[3] <- clx(fm = QuadHolderAidValue,  cluster = outcomesData_long$cas_code)[2, 2]
p_valuesCR_otherAidValue[3] <- clx(fm = QuadHolderAidValue,  cluster = outcomesData_long$cas_code)[2, 4]




##### 	2.5.1 quadratic analysis with controls -- Additional aid value


QuadCovHolderAidValue <- lm(AidValue  ~ outcomesData_long[, "Treatment"] + outcomesData_long[, "WINTELEV_oct13_rescaled"] +  I(outcomesData_long[, "WINTELEV_oct13_rescaled"] * outcomesData_long[, "Treatment"]) + outcomesData_long[, "WINTELEV_oct13_rescaled"]^2  +  I(outcomesData_long[, "Treatment"] * outcomesData_long[, "WINTELEV_oct13_rescaled"]^2) +
                              outcomesData_long$age0to4 +
                              outcomesData_long$age5to12 +
                              outcomesData_long$age13to17 +
                              outcomesData_long$age18to59 +
                              outcomesData_long$age60plus +
                              outcomesData_long$Oct2013Pop +
                              outcomesData_long$disabled +
                              outcomesData_long$origin_gov_0 +
                              outcomesData_long$origin_gov_4 +
                              outcomesData_long$origin_gov_6 +
                              outcomesData_long$origin_gov_10 +
                              outcomesData_long$origin_gov_11 +
                              outcomesData_long$educ_le_priinc +
                              outcomesData_long$educ_pricom +
                              outcomesData_long$educ_middle +
                              outcomesData_long$educ_sec_tech_uni +
                              outcomesData_long$age_head +
                              outcomesData_long$monthinleb +
                              outcomesData_long$base_relatives +
                              outcomesData_long$base_friends 
)
#print(colnames(outcomesData_long)[i])
#print(clx(fm = QuadCovHolderAidValue,  cluster = outcomesData_long$cas_code))
clx(fm = QuadCovHolderAidValue,  cluster = outcomesData_long$cas_code)[2, 1]
clx(fm = QuadCovHolderAidValue,  cluster = outcomesData_long$cas_code)[2, 2]


clx(fm = QuadCovHolderAidValue,  cluster = outcomesData_long$cas_code)

coefs_otherAidValue[4] <- clx(fm = QuadCovHolderAidValue,  cluster = outcomesData_long$cas_code)[2, 1]
sesCR_otherAidValue[4] <- clx(fm = QuadCovHolderAidValue,  cluster = outcomesData_long$cas_code)[2, 2]
p_valuesCR_otherAidValue[4] <- clx(fm = QuadCovHolderAidValue,  cluster = outcomesData_long$cas_code)[2, 4]









# 3. Output Aid Contamination Results Table ----------------------------------


numVar_otherAid <- 2


#  DM inclusion of index and variable selection
RowsNameswSEs_otherAid <- c("$\\hat{\\beta}$: Received Other Aid (yes=1, no=0)", "", "\\emph{p}-value",  
                            "$\\hat{\\beta}$: Total Value of Other Aid (USD)", "", "\\emph{p}-value")



# Create empty vector for control-group means
controlMeans_otherAid <- rep("", times = length(RowsNameswSEs_otherAid))

# Add the control-group means to every 3rd cell
controlMeans_otherAid[seq(from = 1, to = (numVar_otherAid*3)-2, by = 3) ] <- c(round(mean(data$Q46[Treatment==0]), digits = 2), round(mean(AidValue[Treatment==0]), digits = 2))

# Define empty results_otherAid matrix for tex output
results_otherAid <- matrix(nrow = numVar_otherAid*3, ncol = 4)

# add coefficients to tex output
results_otherAid[1, ] <- round(coefs_otherAid, digits = 2)
results_otherAid[4, ] <- round(coefs_otherAidValue, digits = 2)

# add CR SEs to text output
results_otherAid[2, ] <- paste0("(", round(sesCR_otherAid, digits = 2), ")")
results_otherAid[5, ] <- paste0("(", round(sesCR_otherAidValue, digits = 2), ")")



# add FPC CR SEs to tex output
results_otherAid[3, ] <- round(p_valuesCR_otherAid, digits = 3)
results_otherAid[6, ] <- round(p_valuesCR_otherAidValue, digits = 3)



# Bind names, results_otherAid, and control-group means
results_otherAidNames <- cbind(RowsNameswSEs_otherAid, results_otherAid, controlMeans_otherAid)


# Add column names

# Add row at the bottom showing whether Controls were included in each model
results_otherAidNames <- rbind(results_otherAidNames, c("Covariates", "No", "Yes", "No", "Yes", ""))

colnames(results_otherAidNames) <- c(" ", "Linear", "Linear w Covariates", "Quadratic", "Quadratic w Covariates", "Control-Group Mean")

# Output tex table 
print(
  xtable(results_otherAidNames),
  only.contents = TRUE,
  include.rownames = F,
  include.colnames = F,
  floating = F,
  sanitize.text.function = identity,
  sanitize.colnames.function = identity,
  hline.after = c(3, 6),
  file = 'output/otheraid_results.tex'
)





####
# 4. Sharing Results ---------------------------------------------------------
#  53_37_results.tex 


coefs_53 <- c()
sesCR_53 <- c()
p_valuesCR_53 <- c()


coefs_37 <- c()
sesCR_37 <- c()
p_valuesCR_37 <- c()


# 5. Additional Aid Binary ---------------------------------------------------


####  5.1 Linear analysis without controls -- Additional aid binary


LinearHolder53_37 <- lm((data$Q53>0)  ~ outcomesData_long[, "Treatment"] + outcomesData_long[, "WINTELEV_oct13_rescaled"] +  I(outcomesData_long[, "WINTELEV_oct13_rescaled"] * outcomesData_long[, "Treatment"])  )
#print(colnames(outcomesData_long)[i])
#print(clx(fm = LinearHolder53_37,  cluster = outcomesData_long$cas_code))
clx(fm = LinearHolder53_37,  cluster = outcomesData_long$cas_code)[2, 1]
clx(fm = LinearHolder53_37,  cluster = outcomesData_long$cas_code)[2, 2]
clx(fm = LinearHolder53_37,  cluster = outcomesData_long$cas_code)[2, 4]

clx(fm = LinearHolder53_37,  cluster = outcomesData_long$cas_code)

coefs_53[1] <- clx(fm = LinearHolder53_37,  cluster = outcomesData_long$cas_code)[2, 1]
sesCR_53[1] <- clx(fm = LinearHolder53_37,  cluster = outcomesData_long$cas_code)[2, 2]
p_valuesCR_53[1] <- clx(fm = LinearHolder53_37,  cluster = outcomesData_long$cas_code)[2, 4]


##### 	5.2  Linear analysis with controls -- Additional aid binary


LinearCovHolder53_37 <- lm((data$Q53>0)  ~ outcomesData_long[, "Treatment"] + outcomesData_long[, "WINTELEV_oct13_rescaled"] +  I(outcomesData_long[, "WINTELEV_oct13_rescaled"] * outcomesData_long[, "Treatment"]) +
                             outcomesData_long$age0to4 +
                             outcomesData_long$age5to12 +
                             outcomesData_long$age13to17 +
                             outcomesData_long$age18to59 +
                             outcomesData_long$age60plus +
                             outcomesData_long$Oct2013Pop +
                             outcomesData_long$disabled +
                             outcomesData_long$origin_gov_0 +
                             outcomesData_long$origin_gov_4 +
                             outcomesData_long$origin_gov_6 +
                             outcomesData_long$origin_gov_10 +
                             outcomesData_long$origin_gov_11 +
                             outcomesData_long$educ_le_priinc +
                             outcomesData_long$educ_pricom +
                             outcomesData_long$educ_middle +
                             outcomesData_long$educ_sec_tech_uni +
                             outcomesData_long$age_head +
                             outcomesData_long$monthinleb +
                             outcomesData_long$base_relatives +
                             outcomesData_long$base_friends 
)
#print(colnames(outcomesData_long)[i])
#print(clx(fm = LinearCovHolder53_37,  cluster = outcomesData_long$cas_code))
clx(fm = LinearCovHolder53_37,  cluster = outcomesData_long$cas_code)[2, 1]
clx(fm = LinearCovHolder53_37,  cluster = outcomesData_long$cas_code)[2, 2]


coefs_53[2] <- clx(fm = LinearCovHolder53_37,  cluster = outcomesData_long$cas_code)[2, 1]
sesCR_53[2] <- clx(fm = LinearCovHolder53_37,  cluster = outcomesData_long$cas_code)[2, 2]
p_valuesCR_53[2] <- clx(fm = LinearCovHolder53_37,  cluster = outcomesData_long$cas_code)[2, 4]



####  5.3 quadratic analysis without controls -- Additional aid binary 



QuadHolder53_37 <- lm((data$Q53>0)  ~ outcomesData_long[, "Treatment"] + outcomesData_long[, "WINTELEV_oct13_rescaled"] +  I(outcomesData_long[, "WINTELEV_oct13_rescaled"] * outcomesData_long[, "Treatment"]) +  outcomesData_long[, "WINTELEV_oct13_rescaled"]^2  +  I(outcomesData_long[, "Treatment"] * outcomesData_long[, "WINTELEV_oct13_rescaled"]^2) )
#print(colnames(outcomesData_long)[i])
#print(clx(fm = QuadHolder53_37,  cluster = outcomesData_long$cas_code))
clx(fm = QuadHolder53_37,  cluster = outcomesData_long$cas_code)[2, 1]
clx(fm = QuadHolder53_37,  cluster = outcomesData_long$cas_code)[2, 2]

clx(fm = QuadHolder53_37,  cluster = outcomesData_long$cas_code)


coefs_53[3] <- clx(fm = QuadHolder53_37,  cluster = outcomesData_long$cas_code)[2, 1]
sesCR_53[3] <- clx(fm = QuadHolder53_37,  cluster = outcomesData_long$cas_code)[2, 2]
p_valuesCR_53[3] <- clx(fm = QuadHolder53_37,  cluster = outcomesData_long$cas_code)[2, 4]





##### 	5.4 quadratic analysis with controls -- Additional aid binary


QuadCovHolder53_37 <- lm((data$Q53>0)  ~ outcomesData_long[, "Treatment"] + outcomesData_long[, "WINTELEV_oct13_rescaled"] +  I(outcomesData_long[, "WINTELEV_oct13_rescaled"] * outcomesData_long[, "Treatment"]) + outcomesData_long[, "WINTELEV_oct13_rescaled"]^2  +  I(outcomesData_long[, "Treatment"] * outcomesData_long[, "WINTELEV_oct13_rescaled"]^2) +
                           outcomesData_long$age0to4 +
                           outcomesData_long$age5to12 +
                           outcomesData_long$age13to17 +
                           outcomesData_long$age18to59 +
                           outcomesData_long$age60plus +
                           outcomesData_long$Oct2013Pop +
                           outcomesData_long$disabled +
                           outcomesData_long$origin_gov_0 +
                           outcomesData_long$origin_gov_4 +
                           outcomesData_long$origin_gov_6 +
                           outcomesData_long$origin_gov_10 +
                           outcomesData_long$origin_gov_11 +
                           outcomesData_long$educ_le_priinc +
                           outcomesData_long$educ_pricom +
                           outcomesData_long$educ_middle +
                           outcomesData_long$educ_sec_tech_uni +
                           outcomesData_long$age_head +
                           outcomesData_long$monthinleb +
                           outcomesData_long$base_relatives +
                           outcomesData_long$base_friends 
)
#print(colnames(outcomesData_long)[i])
#print(clx(fm = QuadCovHolder53_37,  cluster = outcomesData_long$cas_code))
clx(fm = QuadCovHolder53_37,  cluster = outcomesData_long$cas_code)[2, 1]
clx(fm = QuadCovHolder53_37,  cluster = outcomesData_long$cas_code)[2, 2]


clx(fm = QuadCovHolder53_37,  cluster = outcomesData_long$cas_code)


coefs_53[4] <- clx(fm = QuadCovHolder53_37,  cluster = outcomesData_long$cas_code)[2, 1]
sesCR_53[4] <- clx(fm = QuadCovHolder53_37,  cluster = outcomesData_long$cas_code)[2, 2]
p_valuesCR_53[4] <- clx(fm = QuadCovHolder53_37,  cluster = outcomesData_long$cas_code)[2, 4]






# 5.5  -- Additional aid value

####  5.5.1 Linear analysis without controls -- Additional aid value


LinearHolder53_37AidValue <- lm(data$Q37  ~ outcomesData_long[, "Treatment"] + outcomesData_long[, "WINTELEV_oct13_rescaled"] +  I(outcomesData_long[, "WINTELEV_oct13_rescaled"] * outcomesData_long[, "Treatment"])  )
#print(colnames(outcomesData_long)[i])
#print(clx(fm = LinearHolder53_37AidValue,  cluster = outcomesData_long$cas_code))
clx(fm = LinearHolder53_37AidValue,  cluster = outcomesData_long$cas_code)[2, 1]
clx(fm = LinearHolder53_37AidValue,  cluster = outcomesData_long$cas_code)[2, 2]

clx(fm = LinearHolder53_37AidValue,  cluster = outcomesData_long$cas_code)

coefs_37[1] <- clx(fm = LinearHolder53_37AidValue,  cluster = outcomesData_long$cas_code)[2, 1]
sesCR_37[1] <- clx(fm = LinearHolder53_37AidValue,  cluster = outcomesData_long$cas_code)[2, 2]
p_valuesCR_37[1] <- clx(fm = LinearHolder53_37AidValue,  cluster = outcomesData_long$cas_code)[2, 4]




##### 	5.5.2  Linear analysis with controls -- Additional aid value


LinearCovHolder53_37AidValue <- lm(data$Q37  ~ outcomesData_long[, "Treatment"] + outcomesData_long[, "WINTELEV_oct13_rescaled"] +  I(outcomesData_long[, "WINTELEV_oct13_rescaled"] * outcomesData_long[, "Treatment"]) +
                                     outcomesData_long$age0to4 +
                                     outcomesData_long$age5to12 +
                                     outcomesData_long$age13to17 +
                                     outcomesData_long$age18to59 +
                                     outcomesData_long$age60plus +
                                     outcomesData_long$Oct2013Pop +
                                     outcomesData_long$disabled +
                                     outcomesData_long$origin_gov_0 +
                                     outcomesData_long$origin_gov_4 +
                                     outcomesData_long$origin_gov_6 +
                                     outcomesData_long$origin_gov_10 +
                                     outcomesData_long$origin_gov_11 +
                                     outcomesData_long$educ_le_priinc +
                                     outcomesData_long$educ_pricom +
                                     outcomesData_long$educ_middle +
                                     outcomesData_long$educ_sec_tech_uni +
                                     outcomesData_long$age_head +
                                     outcomesData_long$monthinleb +
                                     outcomesData_long$base_relatives +
                                     outcomesData_long$base_friends 
)
#print(colnames(outcomesData_long)[i])
#print(clx(fm = LinearCovHolder53_37AidValue,  cluster = outcomesData_long$cas_code))
clx(fm = LinearCovHolder53_37AidValue,  cluster = outcomesData_long$cas_code)[2, 1]
clx(fm = LinearCovHolder53_37AidValue,  cluster = outcomesData_long$cas_code)[2, 2]

clx(fm = LinearCovHolder53_37AidValue,  cluster = outcomesData_long$cas_code)

coefs_37[2] <- clx(fm = LinearCovHolder53_37AidValue,  cluster = outcomesData_long$cas_code)[2, 1]
sesCR_37[2] <- clx(fm = LinearCovHolder53_37AidValue,  cluster = outcomesData_long$cas_code)[2, 2]
p_valuesCR_37[2] <- clx(fm = LinearCovHolder53_37AidValue,  cluster = outcomesData_long$cas_code)[2, 4]



####  5.6 quadratic analysis without controls -- Additional aid value 


QuadHolder53_37AidValue <- lm(data$Q37  ~ outcomesData_long[, "Treatment"] + outcomesData_long[, "WINTELEV_oct13_rescaled"] +  I(outcomesData_long[, "WINTELEV_oct13_rescaled"] * outcomesData_long[, "Treatment"]) +  outcomesData_long[, "WINTELEV_oct13_rescaled"]^2  +  I(outcomesData_long[, "Treatment"] * outcomesData_long[, "WINTELEV_oct13_rescaled"]^2) )
#print(colnames(outcomesData_long)[i])
#print(clx(fm = QuadHolder53_37AidValue,  cluster = outcomesData_long$cas_code))
clx(fm = QuadHolder53_37AidValue,  cluster = outcomesData_long$cas_code)[2, 1]
clx(fm = QuadHolder53_37AidValue,  cluster = outcomesData_long$cas_code)[2, 2]

clx(fm = QuadHolder53_37AidValue,  cluster = outcomesData_long$cas_code)



coefs_37[3] <- clx(fm = QuadHolder53_37AidValue,  cluster = outcomesData_long$cas_code)[2, 1]
sesCR_37[3] <- clx(fm = QuadHolder53_37AidValue,  cluster = outcomesData_long$cas_code)[2, 2]
p_valuesCR_37[3] <- clx(fm = QuadHolder53_37AidValue,  cluster = outcomesData_long$cas_code)[2, 4]




##### 	5.6.1 quadratic analysis with controls -- Additional aid value


QuadCovHolder53_37AidValue <- lm(data$Q37  ~ outcomesData_long[, "Treatment"] + outcomesData_long[, "WINTELEV_oct13_rescaled"] +  I(outcomesData_long[, "WINTELEV_oct13_rescaled"] * outcomesData_long[, "Treatment"]) + outcomesData_long[, "WINTELEV_oct13_rescaled"]^2  +  I(outcomesData_long[, "Treatment"] * outcomesData_long[, "WINTELEV_oct13_rescaled"]^2) +
                                   outcomesData_long$age0to4 +
                                   outcomesData_long$age5to12 +
                                   outcomesData_long$age13to17 +
                                   outcomesData_long$age18to59 +
                                   outcomesData_long$age60plus +
                                   outcomesData_long$Oct2013Pop +
                                   outcomesData_long$disabled +
                                   outcomesData_long$origin_gov_0 +
                                   outcomesData_long$origin_gov_4 +
                                   outcomesData_long$origin_gov_6 +
                                   outcomesData_long$origin_gov_10 +
                                   outcomesData_long$origin_gov_11 +
                                   outcomesData_long$educ_le_priinc +
                                   outcomesData_long$educ_pricom +
                                   outcomesData_long$educ_middle +
                                   outcomesData_long$educ_sec_tech_uni +
                                   outcomesData_long$age_head +
                                   outcomesData_long$monthinleb +
                                   outcomesData_long$base_relatives +
                                   outcomesData_long$base_friends 
)
#print(colnames(outcomesData_long)[i])
#print(clx(fm = QuadCovHolder53_37AidValue,  cluster = outcomesData_long$cas_code))
clx(fm = QuadCovHolder53_37AidValue,  cluster = outcomesData_long$cas_code)[2, 1]
clx(fm = QuadCovHolder53_37AidValue,  cluster = outcomesData_long$cas_code)[2, 2]


clx(fm = QuadCovHolder53_37AidValue,  cluster = outcomesData_long$cas_code)

coefs_37[4] <- clx(fm = QuadCovHolder53_37AidValue,  cluster = outcomesData_long$cas_code)[2, 1]
sesCR_37[4] <- clx(fm = QuadCovHolder53_37AidValue,  cluster = outcomesData_long$cas_code)[2, 2]
p_valuesCR_37[4] <- clx(fm = QuadCovHolder53_37AidValue,  cluster = outcomesData_long$cas_code)[2, 4]








# 6. Output Aid Contamination Results Table ----------------------------------


numVar_53_37 <- 2


#  Inclusion of index and variable selection
RowsNameswSEs_53 <- c("$\\hat{\\beta}$: Gave Money to Non-HH Members (yes=1, no=0)", "", "\\emph{p}-value",  
                      "$\\hat{\\beta}$: Helped Other Syrians (yes=1, no=0)", "", "\\emph{p}-value")



# Create empty vector for control-group means
controlMeans_53_37 <- rep("", times = length(RowsNameswSEs_53))

# Add the control-group means to every 3rd cell
controlMeans_53_37[seq(from = 1, to = (numVar_53_37*3)-2, by = 3) ] <- c(round(mean((data$Q53[Treatment==0])>0), digits = 2), round(mean(data$Q37[Treatment==0]), digits = 2))


# Define empty results_53 matrix for tex output
results_53 <- matrix(nrow = numVar_53_37*3, ncol = 4)

# add coefficients to tex output
results_53[1, ] <- round(coefs_53, digits = 2)
results_53[4, ] <- round(coefs_37, digits = 2)

# add CR SEs to text output
results_53[2, ] <- paste0("(", round(sesCR_53, digits = 2), ")")
results_53[5, ] <- paste0("(", round(sesCR_37, digits = 2), ")")



# add FPC CR SEs to tex output
results_53[3, ] <- round(p_valuesCR_53, digits = 3)
results_53[6, ] <- round(p_valuesCR_37, digits = 3)



# Bind names, results_53, and control-group means
results_53Names <- cbind(RowsNameswSEs_53, results_53, controlMeans_53_37)


# Add column names

# Add row at the bottom showing whether Controls were included in each model
results_53Names <- rbind(results_53Names, c("Covariates", "No", "Yes", "No", "Yes", ""))

colnames(results_53Names) <- c(" ", "Linear", "Linear w Covariates", "Quadratic", "Quadratic w Covariates", "Control-Group Mean")

# Output tex table 
print(
  xtable(results_53Names),
  only.contents = TRUE,
  include.rownames = F,
  include.colnames = F,
  floating = F,
  sanitize.text.function = identity,
  sanitize.colnames.function = identity,
  hline.after = c(3, 6),
  file = 'output/53_37_results.tex'
)



# 7. McCrary Tests Section ---------------------------------------------------

### Unadjusted Code
DCdensity(data[, WINTELEV_oct13_rescaled], cutpoint = 0, bin = NULL, bw = NULL, verbose = FALSE, plot = TRUE, ext.out = FALSE)

McCraryPsOct13 <- c()

for(i in   -41:37){

  McCraryPsOct13[i+42] <- DCdensity(data[, WINTELEV_oct13_rescaled], cutpoint = i, bin = NULL, bw = NULL, verbose = FALSE, plot = TRUE, ext.out = FALSE)
  
}

McCraryPMat <- cbind(-41:37, McCraryPsOct13, McCraryPsOct13<=0.05)
# The NAs are the most highly significant differences, which were so close to zero in p-value that they couldn't compute
McCraryPMat[is.na(McCraryPMat[, 3]), 3] <- 1
# Calculate the share of pseudo-thresholds that have significant density discontinuities. 
mean(McCraryPMat[, 3])   

DCdensity(data[, WINTELEV_oct13_rescaled], cutpoint = -41, bin = NULL, bw = NULL, verbose = FALSE, plot = TRUE, ext.out = FALSE)

sum(McCraryPsOct13 <=0.05, na.rm=T)/length(McCraryPsOct13)


McCraryPsJune13 <- c()







#
#
#
#
#
#
#
#
#
###  Code with full registration data base
# Code adjusted alt variable

# Manually identified the minimum altitude value where the function runs
minVal <- -463
#Manually the maximum alt where the function runs
maxVal <- 1471

# range of tests
testRange <- minVal:maxVal

HHAidEligibility[, "WINTELEV_adj"] <- HHAidEligibility[, "WINTELEV"]-500

DCdensity(unlist(HHAidEligibility[, "WINTELEV_adj"]), cutpoint = 0, bin = NULL, bw = NULL, verbose = FALSE, plot = TRUE, ext.out = FALSE)

McCraryPsFull <- c()

HHAidEligibility[, 'WINTELEV_adj' := lapply(.SD, as.numeric), .SDcols = 'WINTELEV_adj'] 

for(i in   1:length(testRange)   ){
  McCraryPsFull[i] <- DCdensity(runvar = unlist(HHAidEligibility[, "WINTELEV_adj"]), cutpoint = testRange[i], bin = NULL, bw = NULL, verbose = FALSE, plot = TRUE, ext.out = FALSE)
}

McCraryPsFull[which(McCraryPsFull=="NaN") ] <- 0

McCraryPMat <- cbind(testRange, McCraryPsFull, McCraryPsFull <=0.05)
# The NAs are the most highly significant differences, which were so close to zero in p-value that they couldn't compute
McCraryPMat[is.na(McCraryPMat[, 3]), 3] <- 0
# Calculate the share of pseudo-thresholds that have significant density discontinuities. 
mean(McCraryPMat[, 3])   


write(round(mean(McCraryPMat[, 3])*100, digits = 2)  ,  file = "output/proxyMcCraryTests.txt")

numberOfRefugeesOct13 <- format((length(unique(HHAidEligibility$ProcessingGroupNumber)) * mean(HHAidEligibility$FS) ), big.mark=",", scientific=FALSE)

write(numberOfRefugeesOct13,  file = "output/numberOfRefugeesOct13.txt")

write( (length(unique(HHAidEligibility$L4DES))  ) ,  file = "output/numberUniqueTown.txt")



write( (length(unique(HHAidEligibility$WINTELEV))  ) ,  file = "output/numberUniqueAlts.txt")



